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The Ginzburg-Landau formalism is constructed for Fermi superfluidity based on the 2n- 
body condensation in parallel with the usual Ginzburg-Landau formalism for the Gooper 
pair condensation. By this formalism, the transition temperature of the 2n-body conden¬ 
sation is given compactly by the zero of coefficient of the quadratic term in the 2n-body 
condensed order parameter without counting complicated Feynman diagrams for the 2n- 
body condensation susceptibility. It is shown that the 2n-body condensed state is stabilized 
in the intermediate- or strong-coupling regime against the BGS state based on the Gooper 
pair condensation. This theory is applied to the case of the possible superfluid state of cold 
atomic gases such as ®Be and ^^^Yb; the quartet state is possible in the former case, while 
the sextet state is expected in the latter case. It is predicted that the ^^^Yb atomic gas so 
far attained satisfies the condition for the sextet condensed state to be realized. 


1. Introduction 

The superfluidity known so far is sustained by the Bose-Einstein or Cooper pair conden¬ 
sate. The former is realized in liquid ^He and some atomic gases of alkali metal elements, 
while the latter is realized in liquid ^He and a variety of superconductors. Over the past 
decade, it has been found that the superfluidity based on the Cooper pair condensation is also 
realized in fermionic atomic gases of alkali metal elements.^) A new aspect of the latter case 
is that the crossover to the Bose-Einstein condensation of diatomic molecules is possible with 
the help of the so-called effect of Eeshbach resonance.®’ 

In principle, there exists another possibility that superfluidity is sustained by a condensate 
based on four fermions (quartet) as in an a-particle correlation in light nucleus.®^ The a- 
particle consists of two protons and two neutrons that have approximate quadruple degeneracy 
corresponding to the 2x2 degeneracy of the real spin and isotopic spin states. In this context, 
the problem of the quartet condensation has been discussed from time to time over the past 
decade or so.®^^®^ The problem of the quartet condensation has also been addressed in the 
context of a fermionic atomic gas with fourfold degeneracy in internal degrees of freedom such 
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as the ®Be atom, which has a nuclear spin I = 3/2 with an electron spin S = 0.^^^ 


The ground state of the four-particle system of such a particle is known to be fully anti¬ 
symmetric with respect to spin coordinates and fully symmetric with respect to space (or 
wavenumber) coordinates.^^) It has been shown, by solving the so-called “Cooper problem”, 
that the quartet state can be stabilized against the Cooper pairing state when four particles 
move outside a rigid Fermi surface in a moderately strong or strong-coupling region of dilute 
systems.It is expected that the quartet superfluid state is possible, in principle, in fermionic 
atomic gases with a nuclear spin 7 = 3/2 and an electron spin S = 0, such as ®Be. 

Such a superfluid state with the 2n-body condensation beyond the Cooper pair condensa¬ 
tion (n = 1) may also be possible for n > 2. The possible n is restricted by the condition that 
2n < 21 + 1 with I being the nuclear spin. For example, the ground state of the ^^^Yb atom 
is sextuply degenerate, i.e., nuclear spin 7 = 5/2 and electron spin 5 = 0, so that a sextet 
condensed state (n = 3) is possible in principle. Indeed, the scattering length analysis (within 
s-wave scattering) of ^^^Yb shows that it is located in a rather strong-coupling region with a 
scattering length a* ~ 11 nm.^®) This implies that the shallow two-body s-wave bound state 
exists, which guarantees the existence of a 6-body bound state because the ground state of a 
6-particle system with 7 = 5/2 is fully symmetric in space coordinates and fully antisymmetric 
in spin coordinates according to the theorem by Nagaoka and Usui,^^) and has a lower energy 
than three 2-body bound states. Namely, in the dilute limit, the 6-body correlation is expected 
to dominate the 2-body correlation, promoting the sextet condensed state of fermionic atomic 
gas of ^^^Yb compared with the Cooper pairing state. However, in the case of an intermediate 
or high density of atoms, these two condensed states compete with each other, as discussed 
in the “Cooper problem” of quartet condensation.^^) Therefore, we need to investigate the 
relative stability of these two states. It was reported that the atomic gas of ^’^^Yb is cooled 
to T/Tp = 0.37 in an optical trap.^^) In these situations, a theory for discussing the 2n-body 
condensation with n > 2 is desired. 

For the quartet condensation, considerable theoretical research studies have been accu¬ 
mulated over the past decade or so, not only as a problem of nuclear physics®^^^) but also as 
a subject of fundamental interest in materials physics. The quartet condensation was shown 
to be possible in one-dimensional models with quadruply degenerate internal degrees of free¬ 
dom. Possible phases of cold atomic systems with a spin 7 = 3/2 were reviewed from 
a wide theoretical point of view.^^) A possibility of four-electron attractive interaction in 
electron-phonon coupled systems was also discussed.^^) However, a concise formalism that en¬ 
ables the estimation of the transition temperature for the 2n-body condensation with n > 2, 
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including the quartet condensation (n = 2), is expected, while some trials have been reported 
in a community of nuclear physics for the quartet condensation®did3) _ 

The purpose of this paper is to construct a Ginzburg-Landau-type formalism for the 2n- 
body condensation with n > 2 in general. It will turn out that this is possible by using 
numerical calculations at a realistic cost for any n as far as the 2n-body condensation with a 
zero center-of-mass momentum is concerned. Namely, a theoretical treatment of d la Nozieres 
and Schmitt-Rink is beyond the scope of the present paper. Nevertheless, we can give a 
physical picture of the 2n-body condensed states. In particular, our result is applicable for 
discussing the possibility of observing the sextet superfluidity in a cold atomic gas of ^^^Yb 
in an optical trap.^^^ 

The organization of the paper is as follows. In Sect. 2, the idea of the Ginzburg-Landau 
(GL) theory for the Gooper pair condensation is extended to the case of the quartet condensa¬ 
tion following the idea of the variational principles of the mean-field approximation based on 
the Feynman inequality for the thermodynamic potential. The GL thermodynamic potential 
(up to quadratic terms) is given in compact form, which is tractable with reasonable computa¬ 
tion time. In Sect. 3, it is shown that the GL formalism is extended to the case of the 2n-body 
condensation without any essential difficulties. Explicit forms of GL thermodynamic potential 
in general form for any n are obtained. In Sect. 4, the transition temperature Tc for the three- 
dimensional free space is calculated for n = 2 ~ 5, namely, from the quartet condensation to 
the dectet condensation, together with the case of the Cooper pair (n = I) condensation. In 
Sect. 5, the case of a two-dimensional square lattice, to which the fast-Fourier-transformation 
(FTT) technique is applicable, is discussed and Tc is calculated for any filling of particles. In 
Sect. 6, the filling dependence of Tc and the quartic terms of the GL thermodynamic potential 
on the square lattice are discussed for the quartet condensation. In Sect. 7, a possibility of 
the sextet condensation in ^^^Yb atomic gas is discussed. In Appendix A, the GL theory is 
reformulated on the basis of the idea of the variational principles of the mean-field approxi¬ 
mation based on the Feynman inequality. In Appendix B, single-particle Green’s function in 
real- and imaginary-time spaces is given. In Appendix C, the expressions of quartic terms in 
GL expansion for the quartet condensation are derived. 

2. Ginzburg-Landau Theory for Quartet Condensation 


In this section, we consider a many-particle system of fermions with fourfold-degenerate 
internal degrees of freedom. For example, ®Be has a nuclear spin 7 = 3/2 and the states with 
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Iz = ±3/2, ±1/3 are degenerate. The Hamiltonian of such a system is expressed as 


®k+q/2,(T®-k+q/2,cr'®-k'+q/2,(T'Ok'+q/2,o-) (1) 

k,(T q k,k' a,a' 

where the summation with respect to the spin variables a and a' is taken over a = 3/2, 
13 = 1 / 2 , 7 = — 1 / 2 , and <5 = —Sj2, which represent the internal degrees of freedom, e.g., Iz- 
Here, £ = k'^/2m being the kinetic energy of particles and fj, being the chemical 

potential, and the two particle interaction Vk^k' is assumed to be independent of the internal 
degrees of freedom. This Hamiltonian is regarded as a generalization of that used in the Cooper 
pair condensation. However, it is more convenient for discussing the quartet condensation or 
the 2n-body condensation with n > 3 to represent Eq. (1) in the form 
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Similarly to the case of the Cooper pair condensation (discussed in Appendix A), the 
mean-field Hamiltonian can be represented as 


Hm{ = Y^'^ 


p^po-^pc 


^ A(pi, P 2 , P3, P 4 ) ^ Pi I al^aaYal^Yp4S + (3) 


P,(T 


Pir - :P4 


vi=l 


where A(pi, p 2 , ps, P 4 ) is the mean field in the present case and assumed to be independent 
of the internal degrees of freedom, a, (3, 7 , and 5, as in the case of the “Cooper problem” 
discussed in Ref. 14, where it was assumed that the wave function for the spin state is fully 
antisymmetric as in the case of the four-particle state. As discussed in Appendix A, the CL 
thermodynamic potential flcL is given explicitly by 


f^GL—± {kl IImf)mf- 

The operator corresponding to the second term in Eq. (4) is given by 
= ^ ^Hp4_P45(pi-f P2-P3-P4)aj,i,a( 


( 4 ) 


P2,^‘ 


'^P 3 ,cr'^k 4 ,cr 


PIT" jPt a,a' 

+ ^ A(pi,P 2,P3,P4)<5 a|,4„aj,^^at3^aj,^5 + h.c.. 


( 5 ) 


Pl) - )P4 


7=1 


First, we calculate Hmf by perturbation expansion with respect to the mean field A in 
the mean-field Hamiltonian (3) up to the quadratic term in A(pi, p 2 , ps, P 4 ) because we are 
interested in obtaining the transition temperature for the moment. The quartic term will be 
discussed later. Hereafter, we assume that the wavenumber dependence of A is fully symmetric 
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with respect to p* (i = 1 ~ 4) as in the case of the four-particle ground state,and is given 
with a variational function /(p) as 

4 

A(pi, P2, P3, P4) = A /(pi). (6) 

i=l 

The wave function /(p) is a generalization of that introduced in Ref. 14 for the “Cooper 
problem” of the quartet bound state. 

The result for is given as 

~flo-A4(r)|A|2 + 0(|A|4), (7) 


where flo is the thermodynamic potential in the normal state, and the coefficient A 4 ^{T) is 
given by the Feynman diagram shown in Fig. 1, and its analytical expression is given as 


y44(r) = n/ ^ J]l/(R)pG'(pi,ie„J 

i=l €ni 

X(5(pi +P2 +P3 +P4) X +, +e +e„,,0, 


( 8 ) 


where G is the Matsubara Green function of quasiparticles in the normal state and is assumed 
to be independent of the four spin variables a, /3, 7 , and S. Hereafter, = (2n + l) 7 rT is the 
fermionic Matsubara frequency. By using the identities 

5 (Pl + P2 + P3 + P4) = I ^e'(Pl+P2+P3+P4)-r^ (9) 

and 


€n H”^n2 


Q = T / (^.,-g-lFni+':n2+en3-|-en4)T^ 

Jo 


the coefficient A^iT)^ given in Eq. ( 8 ), is reduced to a compact form as 




dp 

W)' 


T |/(p)pG(p,ien)e'(P' 


r-e„r) 


( 10 ) 


( 11 ) 


Note here that the numbers of integration and summation variables are greatly reduced. This 
point is much more crucial for extending the discussion to the cases of the sextet, octet, and 
dectet condensations. 

Next, we calculate the grand canonical average of Eq. (5) with the mean-field Hamiltonian 
(3) up to quadratic terms in the gap A. These terms are given by the Eeynman diagrams shown 
in two terms of Eig. 1 (with a positive sign) and Eig. 2, and their analytical expressions are 
given as 


{H - - 2yl4(T)|A|2 + EH4(r)|A|2 + 0(|A|4), (12) 


where we have assumed that the two-particle interaction is wave-vector-independent with 
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Pi, a 



A* 


Fig. 1. Feynman diagram for f2„if of the quadratic term with respect to A and A*. 


the energy cutoff Eq (on the order of the Fermi energy ep) considering the case of a dilute 
atomic gas with an s-wave attractive interaction, or a model case of fermion with multi¬ 
internal degrees of freedom moving on a lattice. The expression for B/i{T) in Eq. (12) is given 
as 

o 

dpi 


B^{T) = iC2T^X[ 


2=1 ' 


(2vr) 
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1=3- 


(27r)' 
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dp' 
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X(5 (pi + P2 + P3 + P4) X 

X(5 (pi + P 2 + p's + P 4 ) X (13) 

Here, the combination factor 46*2 represents the number of ways of choosing two (connected 
to the interaction V) of four Green functions. By using Eqs. (9) and (10) and similar ones, 
the coefficient B 4 {T) is reduced to 

2 ,/3 


H 2 


B. 


2 — 1 

/ (^^^IZ-^(‘l2)G'(q2,iem2)e'(‘i2' 


ri) 


■r2—em2T2) 


2 


dq3 

(27r)’ 


Tj]|/(q3)pG(q3,ie,n3: 


oi(q 3 '(ri+r 2 ) £7713 (ri-|-T2)) 


(14) 


This expression is also numerically tractable as that for A 2 n{T), given in Eq. (11). This is also 
the case for the sextet, octet and dectet condensations, as discussed in the next section. 
Adding Eqs. (7) and (12), the GL thermodynamic potential Hgl(^) is expressed as 

Hgl(A) ~ Ho + [A4(r) + VB4{T)] jAp + 0{\A\^). (15) 
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Pi, a 


A 



P4,<5 V P4,<5 


A* 


Fig. 2. Feynman diagram for {H — including the interaction V of the quadratic term with 


respect to A and A* . 


Then, the transition temperature for the quartet condensation is determined by the relation 


Ai{T,) = \V\B^{T,). 


(16) 


This is a natural extension of that for the Cooper pair condensation, i.e., Eq. (A-12), leading 
to the BCS formula Eq. (A-13). 

The quartic terms in A and A* include the integration with respect to ri, r 2 , and ra, 
and Ti, r 2 , and rs, respectively, as discussed in Appendix C [See, e.g., Eqs. (C-4) and (C-9)]. 
Therefore, integrations similar to Eqs. (11) and (14) are technically impossible to perform 
within a reasonable computation time in the case of three-dimensional space, which will be 
discussed in Sect.4. 

On the other hand, in the case of a two-dimensional square lattice, it is possible to perform 
the calculations by exploiting the technique of fast Eourier transformation (EET), as will 
be discussed in Sect.6, in which the quartic term will be shown to be positive for relevant 
parameter sets. Therefore, it is reasonable to assume that the quartic term with respect to A 
and A* has a positive finite value also in the case of three-dimensional free space , making 
the transition a second-order one. 

3. Generalization to 2n-Body Condensation 

The formalism determining the transition temperature developed in the previous section 
for the quartet condensation is easily generalized to the sextet, octet, and dectet condensations. 

Eor the sextet condensation, the coefficient Aq{T) is given by the Eeynman diagram shown 
in Eig. 3, and its analytical expression is given in parallel with Eq. (8) as follows: 



x5 (pi -b P2 + P3 + P4 + P5 + Pe) X (17) 


where G is the Matsubara Green function of quasiparticles in the normal state and assumed 
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to be independent of the six spin variables a, 13, 7 , C) 7 , and By using the identities, similar 
to Eqs. (9) and (10), 


<5 (Pl + P2 + P3 + P4 + P5 + P6) = j ^eHPl+P2+P3+P4+P5+P6) r^ (Ig) 

and 


“rtn2 “rtrig 


Q = X —i(*:rij+en2+CTi3+en4+en5+eng)'r 

'’Jo 


the coefficient Aq{T), given in Eq. (17), is again reduced to a compact form as 


A, 


r I i/(p)po(p. 


) = 


r-enT) 


(19) 


( 20 ) 


The numerical calculation of Eq. (20) can be performed at the same computational cost as 
Eq. (11). Namely, the increase in the integral or summation variables in Eq. (17), compared 
with that in the case of quartet condensation, is absorbed by the identities Eqs. (18) and (19). 


Pi,a 



A* 


Fig. 3. Feynman diagram for Bmf of the quadratic term with respect to A and A*. 


Similarly, the coefficient Bq{T), whose Feynman diagram is given by Fig. 4, is calculated 
in parallel with Eq. (13) as follows: 


B,iT) = oC2T^fl I ^ J]|/(p*)pG(p„i6„J 


6 


"'ll / J^[/(Pi)]V(p;)G'(p,-,ie„,)G(p',ie, 

i=5"' '' ’ e„. ^ ‘ ■ 
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(27rl3 

X(^ (P1 + P2 + P3 + P4 + P5 + Pe) X 

X J (pi + P2 + P3 + P4 + P 5 + Pe) X ,0. (21) 

Here, the combination factor qC 2 represents the number of ways fo choosing two (connected 


to the interaction V) of six Green functions. By using Eqs. (18) and (19) and similar ones. 
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2 



2 


2 




4 



( 22 ) 


The calculation of Eq. (22) is performed at the same numerical cost as Eq. (14) for the quartet 
condensation. 

Then, the transition temperature of the sextet condensation is also given by Eq. (16) 
with Aq{T), given by Eq. (20), and Bq{T), given by Eq. (22), as in the case of the quartet 
condensation. 


Pi,a 



A 


A* 


Pfi,? V 


Fig. 4. Feynman diagram for {H — including the interaction V of the quadratic term with 


respect to A and A*. 

As one can see from the derivation of the coefficients Aq{T) and Bq{T) above, one can infer 
a general expression for these coefficients. Namely, for the octet condensation, the exponent 
of [J (ip/(27r)^ T •••] in the expression of Aq{T), i.e., Eq. (20), is only replaced by 8, 
and the exponent of the last factor [J dq3/(27r)^ T •••] in the expression of Be(T), 

i.e., Eq. (22), is only replaced by 6. This is easily generalized to the case of higher number 
of condensation unit, say the octet or dectet condensation. Eor 2n-body condensation, the 
exponent of [/ (ip/(27r)^ ''' ] ™ expression of A 2 n(T), i.e., Eq. (11), is given by 2n, 

and the exponent of the last factor [J fiq3/(27r)^ • • • ] in th® expression of B 2 n(T), i.e., 

Eq. (14), is given by 2(n — 1). A combination factor of B 2 n{T) is given by the number of ways 
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of choosing 2 lines from 2n lines, i.e., 2 nC '2 in general. Namely, A 2 n{T) and B 2 n{T) are given 
by the following expressions: 

P \ -["in 

A2n{T) = J dT y’^^r^|/(p)|2G(p,ie„)e'(P~") , (23) 

L 

and 

2 r 1 ^ 

^2n(r) =dn J dri J 

2=1 €777,^ 

1 ^ 

€7712 

-1 2(n-l) 

X J -( 24 ) 

emg 

Then, the transition temperature of 2n-body condensation is determined by Eq. (16) by 
using the coefficients A 2 n(T) and i? 2 n(^) instead of vl 4 (T) and B^iT). It is remarkable that 
the numerical cost for the transition temperature Tc does not increase with increasing number 
n for 2n-body condensation. This is a secret of attacking the problem using the generalized 
Ginzburg-Landau formalism. 


4. Three-Dimensional Free Space 

In this section, we calculate the transition temperatures for 2n-body (n = 2 ~ 5) con¬ 
densation and compare them with that for the Cooper pair condensation in three-dimensional 
free space. Precisely speaking, the variational wave function /(p) in Eq. (6) should be deter¬ 
mined so as to minimize the thermodynamic potential or free energy. However, since such a 
calculation needs a much longer time, we here adopt an approximate solution by assuming 


/(P) 


/ 

I , 0 ^ Ep ^ Eq , 

< 

0, Eq <i Ep. 


(25) 


Nevertheless, a fundamental aspect of 2n-body condensation is expected to be captured. 

Let us define the quantity in the square brackets of Eq. (II) by G'(r, r) which is given 
explicitly as follows (see Appendix B for its derivation): 

m ._ p{p-r)£F^ 

T y_i ^^sin[^/^(/cFr)] ( 26 ) 

where ep and fcp are the Fermi energy and Fermi wave number, respectively, and Xq = {eq/e-p) — 
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1. The x-integration in Eq. (26) well converges as Xc increases. Then, we put Ec = 2eF, i.e., 
Xc = 1. It turns out by explicit numerical calculations that the integrations with respect to r in 
the expressions of A 2 n{T), i.e., Eqs. (11) and (20), and B 2 n{T), i.e., Eqs. (14) and (22), should 
be taken over a sufficiently wide r-region. On the other hand, angular integration with respect 
to the direction r = r/r is easily performed, giving only the factor dvr. Therefore, proper r- 
integration remains to be performed. In order for A 2 {T) for the Cooper pair condensation to 
exhibit a logarithmic T dependence down to T = 10“^ epj we have to take the integration 
over 0 < r/cp < 1000/\/2. Moreover, the contribution from the region rk-p <C 1 should also 
be calculated properly so that we have to take finer meshes there. Therefore, we choose the 
following points on the r-axis: 


- I 

rn=r*P^, (27) 

Pr- I 

and take the summation from n = 1 to n = by multiplying the width of each mesh, 
Ari = r*(l + pr/‘i) for n = 1, and 


= 


rn + Vn+l rn-l+Vn r 

-^-^= Ur + l), 


(28) 


for 2 < n < Mr- Namely, we use a modified trapezoidal rule. Explicitly, we take r* = 
10“®/(\/2A:f), Pr = 1.02419764544894, and Mr = 1000, which yields rMrt given by Eq. (27), 
VMr - 1000.0000000/(\/2A:f). 


On the other hand, the r dependence of G{y, r), given by Eq. (26), near r = 0 and j3 is very 
sharp in the limit fiep S> 1 because the factor _j_ exponentially small for 

— 1 < X < Xc in the intermediate region 0 < r < /3, while it is nearly equal to 1 for 0 < x < Xc 
and —l<x<0atr = 0 and r = /3, respectively. Therefore, it is crucial to properly take into 
account the sharp variation of G{t, t) near r = 0 and /3 in numerical integrations in Eqs. (11), 
(14), (20), and (22). To this end, we take meshes of the r-integration as follows. Similarly to 
the case of r-integration, we choose the following points in 0 < r < /3/2 on the r-axis 

./>r-l 


Tn = n + T 


and in /3/2 <T<j3 


Tn= P-TS-T 


Pr-l' 
Pr-^ 


(29) 


(30) 


Pr-l 

and take the summation from n = 0 to n = [M^- being chosen as an even natural 

integer) by multiplying the width of each mesh: Atq = + t* j2 for n = 0, and 


At„ = + 1 ). 


(31) 
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for 1 < n < Mr/2, and 

= 'r(Mr/2) 


-TjMrm-l + T-^Mr/2) ^ ^(Mr/2)-l 

2 2 ’ 


(32) 


for n = Mr/2. Here, we have introduced a small ts in order to avoid singular behaviors at 
r = 0 and t = j3. Explicitly, we take t* = 10“®/3, ts = 10“®/?, Pr = 1.49932125806831, and 
Mr = 50, which yields T(m,-/ 2 )) given by Eqs. (29) and (30), Ti^Mrli) — 0.5000000/?. 

The relation determining the transition temperature T^, i.e., Eq. (16), is transformed to 


1 = |i^lx2„(rc), 


(33) 


where the “2n-body condensation susceptibility” X 2 n{T) is defined by 

_ B2n{T) 

A2n{Ty ^ ^ 

where A 2 n{T) and B 2 n{T) are the expressions for 2n-body condensation, respectively: e.g., 
A 2 n{T) and i? 2 „(T) are given by Eqs. (11) and (20), and Eqs. (14) and (22) in the cases of the 
quartet (n = 2) and sextet (n = 3) condensations, respectively. However, it should be noted 
that X 2 n{T) cannot be represented by a canonical correlation function of any quantities. This 
is in marked contrast with the case of the Cooper pair condensation (n = 1), in which X 2 {T) 
is given by A 2 {T) whose explicit form is given by 


A2{T) = 


dr 

(M 


dr 


dp 


n 2 


T ^ |/(p)pG(p,ie„ 


„i(p.r-e„r) 


(35) 


with the same energy cutoff £c = as that in the case of n > 2. This A 2 {T) is simply Ki{T), 
given by Eq. (A-7), which is the canonical correlation of the pair operator, as discussed in 
Appendix A. 

Eigure 5 shows the temperature dependence of the “2n-body condensation susceptibility” 
X2n{T) for n = 2 ~ 5, i.e., from the quartet condensation to the dectet condensation, together 
with the Cooper pair susceptibility X 2 {T) = A 2 {T) =Ki{T). The unit of X 2 n is N{£-p), the 
density of states at the Eermi level per spin component. This result implies that the 2n-body 
condensation (with n > 2) has a larger “susceptibility” than the Cooper pair condensation 
in the high-temperature region T > 10“^eF and vice versa. Another intriguing aspect is that 
there exists a threshold coupling, |hth|, necessary for 2n-body condensation to occur at T = 
0 K, and a reentrant of the superfluid state is expected as the temperature T is decreased in 
the case of |E| > |V/;h|. On the other hand, if E < 0, the Cooper pair condensation is always 
possible as T is sufficiently reduced, no matter how the is low, because X 2 {T) = Ki{T) 
diverges logarithmically in the limit T —)■ 0. This is consistent with the result for the stability 
of the quartet condensation against the Cooper pair condensation at the level of the “Cooper 
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problem” discussed in Ref. 14 , in which the quartet state has a lower energy than two Cooper 
pairs only in the intermediate- or strong-coupling region. 



Fig. 5. (Color) “2n-body condensation susceptibility” X 2 n{T)/N{e-p) as a function of temperature 
T/ep in a logarithmic scale. 


Here, we discuss why X 2 n{T)'s (n > 2) exhibit peaks at T ~ (0.2 ~ 0.4)eF, as shown in Fig. 
5. For an explicit discussion, we discuss the case of the quartet (n = 2) condensation. First, we 
note that H 4 (T), given by Eq. (8), is the “bare” susceptibility of the quartet condensation, as 
shown in Fig. 1. The T dependence of A^^T) is shown in Fig. 6, in which one can see that Ai[T) 
exhibits a peak at T ~ 0.2eF. Therefore, the quartet susceptibility has a tendency of exhibiting 
a peak structure at around T = 0.2eF. In the high-T region, T 3> ef, A 4 {T) oc T~^, so that 
Aji{T) increases as T decreases T = e-p because a restriction on momentum integrations due 
to the momentum conservation law in Eq. (8) is less severe in the classical region (T > sf) 
than in the Fermi degenerate region (T <C sf). On the other hand, in the low-T region, 
i.e., T <C ef, A 4 ^{T) should decrease (to a certain finite value) as T decreases because the 
restriction due to the momentum conservation law becomes crucial owing to the effect of 
Fermi degeneracy, which suppresses the available momentum space. As a result, the peak 
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structure in Ai{T) is expected to appear. This is in marked contrast to the case of the Cooper 
pair condensation, for which X 2 {T) is given by A 2 {T), given by Eq. (35). Since A 2 {T) is 
free from such an extra restriction due to the momentum conservation law, A 2 {T) increases 
monotonically (logarithmically) as T decreases . Similarly, given by Eq. (13), appearing 

in the numerator of X 4 (E), given by Eq. (34), also exhibits a more pronounced peak structure 
than A4^{T), as shown in Eig. 6 . This is because oc at T ep so that i?4(T') 

increases more sharply than ^ 4 (T) as T decreases , making the peak height much higher. As 
a result, a peak structure in = Bi{T)/A 4 ^{T) appears at around T = 0.2eF. 



Fig. 6. (Color online) A 4 ^{T)/N{e-p), B 4 ^{T)/[N{ep)Y, and Xi[T)/N{ep) as a function of temperature 
T/ep in a logarithmic scale. 


The transition temperature Tc/ep determined by Eq. (33) is shown in Eig. 7 as a function 
of the strength of attractive interaction |E|Ai(ep) for a series of 2n-body condensations. The 
Tc of 2n-body condensation (n > 2) is higher than that of the Cooper pair condensation in 
the intermediate-coupling region, |E|A^(ep) < 1 and strong-coupling region, |E|A^(ep) > 1. 
Namely, the 2n-body (n > 2) condensed state is stabilized against the Cooper pair condensed 
state in such regions. Eor the attractive interaction V ~ Ej;h) Ec exhibits a reentrant behavior. 
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However, such a region of |H| is restricted in a very narrow region above |Vth|- The threshold 
strengths of V are |Tth|-/V(eF) — 0-1, 0.15, 0.2, and 0.3 for the dectet, octet, sextet, and 
quartet condensations, respectively. Indeed, in a wide region \V\ > |Hth|) 2n-body (n > 2) 
condensations dominate the Cooper pair condensation. 

On the other hand, in the strong-coupling region |H|A^(eF) > 1, we have to take into 
account the effect of the center-of-mass motion of such molecules beyond the mean field 
approximation adopted in previous sections, in which the center of mass is assumed to be 
at rest. Then, Tc is determined by the Bose-Einstein condensation temperature TbeC; which 
is higher in the case of diatomic molecules than in the case of 2n-atomic molecules. This is 
because the mass of a 2n-atomic molecules is n times larger than that of a diatomic molecule, 
and the number density N/2n of 2n-atomic molecules is 1/n times smaller than that of 
diatomic molecules, resulting in the Tc of a 2n-atomic molecule gas being times smaller 

than that of a diatomic molecule gas, since Tbec is given as A:bTbec (/iV^b) X (iVb/E)2/3, 
rub and N\y/V being the mass and number density of a composite boson. In this strong¬ 
coupling region, we need to extend the theory so as to take into account the center-of-mass 
motion, as in the theory of Nozieres and Schmitt-Rink for the BCS-BEC crossover of the 
transition temperature.However, this is beyond the scope of the present study, and is left 
for future studies. 

5. Two-Dimensional Square Lattice 

In this section, we discuss the problem in the two-dimensional tight binding model on 
the square lattice with nearest-neighbor transfer. The energy dispersion of this model is well 
known: 

£k = —2t(cos kxa -|- cos kyo), (36) 

where t is the transfer integral among nearest-neighbor sites and a is the lattice constant. In 
the lattice model, the attractive interaction at the on-site is denoted as —[/, which should 
be distinguished from the Fourier component V of the interaction in the continuum model in 
three-dimensional free space discussed in previous sections and Appendix A. Corresponding 
to Eq. (26), the Matsubara Green function G(rj,r) at the lattice point r* and the imaginary 
time r is given by 

k en 

where Nj, is the number of lattice points and — /x. Note that G(rj, r) is a real quantity 

because it is given by an expression similar to Eq. (B-4), which is real since the term including 
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Fig. 7. (Color) Phase diagram in the Tc/ef — |F|A^(eF) plane for the 2n-body condensation (n > 2) 
and the Cooper pair condensation. 


cos(p • Tj) vanishes. In order to apply the technique of fast Fourier transformation (FFT) to 
the calculation of the coefficients A 2 n{T) and B 2 n{T) given in Sect. 2, let us introduce the 
following quantity: 

Xm{ri,T) = [G{ri,T)Y^. (38) 

Note that Xm(rj, r) is a real quantity and expressed by the Fourier series as (in the case where 
m is an even natural number) 

Xm{v^,r) = ^ (39) 

where the Fourier component Xm(kj,icjn) is defined as 

Xm(kj,ia;n) = ^ / dr^^(ri,(40) 

^ Jo 

where ujn = 2TmT is the bosonic Matsubara frequency because Xm{'ri,T + (3) = Xm(rj,r). 
The coefficients A 2 n{T), given by Eqs. (11) and (20), and B 2 n{T),given by Eqs. (14) and (22), 
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are expressed in terms of Xm(ri,r), given by Eq. (38), as follows: 


A2n{T) = E/ (irX2„(ri,r), 


(41) 


and 

B 


2n{T) = 2nC2 EE/ dn f (ir2X2(rf\ri)X2(rf\r2)X2n-2(rf^+rf\ri+ r2). ( 42 ) 

nt (o\ */0 


r(l) 

^ i ^ i 


Substituting Eq. (39) into Eqs. (41) and (42), and taking summations with respect to r*, 
and and performing integration with respect to r, ri, and r 2 , these quantities are 
expressed in terms of the Eourier component in Eq. (40) as 


A 2 n{T) — ^ 2 n (kj, icUri 


(43) 


k-,’ 


and 


B2n{T) — 2nC2 ^ ^ [X2(—kj, —icu^)]^ X2„_2(kj, iWn). (44) 

^ kj UJn 

A number of k-points in the two-dimensional Brillouin zone is taken as 2^ x 2^, and that 
of the bosonic (fermionic) Matsubara frequency Un = 2mrT (cn = (2n -|- l)vrT) is restricted 
within the region —2^® < u < 2^®. One may suspect that this number of meshes 2^ x 2^ is 
not sufficiently large to maintain the accuracy of the results. However, we have verified that 
this number gives sufficient accuracy by performing calculations for a series of numbers of 
meshes by relaxing the cut in the Matsubara frequency, which is much more important for 
maintaining the accuracy of calculations. Nevertheless, this mesh size gives a restriction on 
temperature above which the accuracy of calculations of A 2 n{T) and H 2 „(T) is guaranteed. 
The lower limit of the temperature Tll is estimated as follows: Tll = 8t/(2^ x 2^), where 8t 
is the bandwidth of dispersion of Eq. (36) and 2^ x 2^ is the number of meshes in the first 
quadrant in the Brillouin zone. This restriction for temperature, T > Tll, is expected to give 
a more severe effect in the case with a low filling of particles compared with half-filling. 

Then, we only have to perform summations with respect to and r or k and Un (or 
€n) several times, instead of directly performing multiple integrations with respect to and 
Tj. The latter calculation needs a much longer time than the present EET technique, and it 
is technically impossible to use it for integrations and summations for Ci{T) [i = 1,2) and 
D^{T) {i = 1~ 9), which are the coefficients of the quartic terms in A and A*, as discussed 
in Appendix C. 

The transition temperature of “2n-body condensation” is given by Eq. (33) with the 
“2n-body condensation susceptibility” X 2 n{T), given by Eq. (34). Eigure 8 shows the tem- 
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perature dependence of X 2 n{T) in the cases from the quartet (n = 2) condensation to the 
dectet (n = 5) condensation together with the case of the Cooper pair condensation (n = 1). 
Owing to a restriction on the size of the number of Matsubara frequencies, there exists a lower 
limit of temperature, Tll, below which the FFT calculation becomes inaccurate. Therefore, 
we show X 2 n{T) for T > Tll in Fig. 8. The filling of fermionic atoms is fixed at the half-filling 
(riA = 1). Here, the filling ua is defined by the ratio of twice the number of occupied states 
in the k-space (in the hypothetical normal ground state) and the total number of k points in 
the Brillouin zone. 

Note that X2n(7")’s for n > 2 have peaks at around T ~ t, and are larger than that 
for n = 1 (Cooper pair susceptibility) in the high-temperature region T > while the 

tendency is reversed in the low-temperature region, i.e., the Cooper pair susceptibility X 2 
dominates X 2 n{T) for n > 2 at T < This is consistent with the result in the case of 

three-dimensional free space shown in Fig. 5. Also note that the combination factor 2 nC 2 is 
crucial for X 2 n{T) with n > 2, which exceeds X 2 {T) in the high-temperature region. Indeed, 
without the factor 2 nC 2 , X 2 {T) is larger than X 2 n{T) with n > 2 in the entire temperature 
region, although we do not explicitly show the result here. 



Fig. 8. (Color) Temperature dependence of “2n-body condensation susceptibility” X 2 n{T) with n>2, 
and the Cooper pair susceptibility at half-filling. The units of energy and temperature are chosen 
as t, the transfer integral among the nearest-neighbor sites. 


Tll is estimated as follows. The maximum magnitude of the Matsubara frequencies is 
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27r X 2^*^ T. Tll is defined by the condition 2ti x 2^° Tll = 40t, where 40t is 10 times half the 
bandwidth 4t, i.e., Tll = 40t/[27r x 2^°] ~ 6.2 x 

Figure 9 shows the relationships between the strength U of the attractive interaction 
and the transition temperature Tc, which is also obtained by solving Eq. (33) in the case 
of half-filling. Here, we show only Tc such that Tc > Tll, as in Fig. 8. In order for the 
“2n-body condensation” with n > 2 to appear, the attractive interaction needs to exceed a 
threshold, while the Cooper pair condensation is always possible, if the temperature is reduced 
sufficiently, owing to a logarithmic divergence of X2(T) oc — log T in the limit T —>• 0. This 
behavior is also consistent with the result in the case of three-dimensional free space shown 
in Fig. 7. 



6. Properties of Quartet Condensation on Square Lattice 

In this section, some aspects of the quartet condensation on the square lattice are dis¬ 
cussed. All the calculations in this section are performed by taking into account the T depen- 


19/48 








J. Phys. Soc. Jpn. 


Full Paper 


dence of the chemical potential /U in a noninteracting system. 

6.1 Dependence on filling of fermionic atom 

Figure 10(a) shows the temperature dependence of Xq{T) [= X^{T)], and Fig. 10(b) shows 
the relationship between the transition temperature Tc and the strength of attractive inter¬ 
action, U/t, for the quartet condensation for a series of fillings ua of fermionic atoms. This 
result implies that Tc increases as the filling increases, which is consistent with the results in 
Refs. 8 and 11. 

Figure 11 shows the relationship between ua and Tc/t of the quartet condensation (shown 
by dots) together with that of the Cooper pair condensation (shown by lines) for a series of 
strengths U of the attractive interaction. One can see that the region with the condensation 
extends to the region of low density (u-a) as U increases. For U/t > 2.5, the Tc of the quartet 
condensation is higher than that of the Cooper pair condensation for any filling 0 < ua < 1- 
This is consistent with the result of the “Cooper problem” in the quartet case, in which the 
quartet state is stabilized in the intermediate- or strong-coupling region and in the low-density 
region,and also consistent with those for the Tc of the a-condensation in the nuclear matter 
discussed in Refs. 8 and 11. On the other hand, in the case of weak and intermediate couplings 
U/t < 2.25, the condensed state appears only in the region ua > n^A^ where denotes a 
threshold filling, and Tc exhibits a reentrant behavior near the threshold ua > • This is 

somewhat different from the results shown in Refs. 8 and 11, where the Tc of the Cooper pair 
condensed state is higher than that of the quartet state in the high-density region, and also 
from the result for the “Cooper problem” discussed in Ref. 14. 


6.2 Quartet ordered state in GL region 

By extending the expression (15) in Sect. 2, the GL free energy of the quartet condensation 
is given in its usual form as 

ii(A) flo + a(T)|A|2 + ]^b{T)\/\\‘^ + • • • , (45) 

where the coefficients a and b are defined as 


and 


a{T) = [A^{T)+V BfiT)] = Ai{T)[l + Vxa{T)] 




VA^{Tc) 


dX4{T) 

dT 


T=Tc 


X {T-Tc), 


(46) 


b{T) = bi{T) + Vb2{T), 


(47) 
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Fig. 10. (Color) (a) T dependence of xq(T), and (b) Tc{Q)/t vs U/t, for a series of fillings riA = 
0.2, 0.4, 0.6, 0.8, 1.0. Data in both figures are restricted to those at T > Tll- 


with bi{T) and h 2 {T) defined as 

hi{T) 

hiT) 


3 


2^[C,(r) + A(T)], 


i=l 

9 


2E 


D,(T) 

V 


(48) 

(49) 
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U/t=3.00 

U/t=2.50 

U/t=2.00 

U/t=1.50 

U/t=1.00 


Fig. 11. (Color) TiA vs Tc for quartet (dots) and Cooper pair (lines) condensations for a series of 
strengths of attractive interaction U/t. Data are restricted to those at T > Tll- 


where Ci{T) (i = 1 ~ 3) and DiiT) (i = 1 ~ 9) are explicitly given in Appendix C. Note that 
the interaction V is equal to —U in Sect. 6.1. It turns out that b(T) is positive by explicit 
calculations below. Therefore, the standard treatment for the second-order phase transition 
is possible. 

Indeed, C'i(T)’s and Zli(r)’s are calculated in Appendix C as follows: C'i(T) and 6 * 2 (T) = 
6 * 3 (r) are given by Eqs. (C- 6 ) and (CTl), respectively; Di{T) = —2Ci{T) and D 2 {T) = 
D^{T) = — 2 C' 2 (r); other coefficients D 4 {T), D^{T), Dq{T), D^iT), Ds{T), and Dg{T) are 
given by Eqs. (C-23), (C-27), (C-32), (C-36), (C-40), and (C-44), respectively. The results of the 
filling (chemical potential /r) dependences of bi{T)t^ and b 2 {T)t‘^ at T = O.lt and T = t are 
shown in Eig. 12. Meshes of summations in these formulas are taken as 2® x 2® for summations 
in wave numbers over the whole Brillouin zone of the square lattice, and as 2® for those in the 
Matsubara frequencies —27t 2‘^T < LUm < 27r2'^T and — 7 r (2 x 2^ — 1)T < < 7 r (2 x 2^ — l)r. 

A lower limit of temperature, above which the accuracy of calculations is guaranteed, is 
dehned by the condition 2Tr2^T^^ = 40t as in Sect. 5, i.e., = 40f/2®7r ~ 4.0x 10“^f. We have 

verified in the case of /r/t = 0 that the accuracy of the temperature dependences of 6 i(T) 
and 62 (E) is maintained up to 90% of those obtained for meshes 2® x 2® for summations in 
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wave numbers and 2® for those in the Matsubara frequencies, which corresponds to the lower 
limit of temperature of = 5.0 x 

One can see in Fig. 12 that b{T) = bi{T) + Vb 2 {T) is positive, at least in the region 
of attractive interaction giving < t (see Fig. 9). Therefore, the phase transition is of the 
second kind, as in the case of the Cooper pair condensation. The results of the temperature 
(T) dependences of bi{T)t^ and b 2 {T)t^ are shown in Fig. 13 for the filling corresponding to 
/x/t = 0 and /x/t = 3.9. This also shows that b{T) = bi{T) + Vb 2 {T) is positive for relevant 
parameter sets giving a reasonable Tc, as shown in Fig. 9, guaranteeing the second-order phase 
transition. 

Of course, there is no technical difficulty in calculating b{Tc) with Tc determined by the 
condition a{Tc) = 0, i.e., ^ 4 (Tc) = \V\B 4 ^{Tc), in Eq. (46). The coefficients bi{T) and b 2 {T) of 
the quartic term in A and A* can be calculated with the required accuracy. Thermodynamic 
analysis based on the GL thermodynamic potential is left for future studies. 



Fig. 


12. (Color online) Filling (chemical potential /x) dependences of the coefficients bi{T)t^ and 
b 2 {T) at T = 0.1 1 and T = t. 


7 . Possibility of Sextet Condensation in Atomic Gas 

It has been reported that ^^^Yb atomic gas is cooled down below the Fermi degeneracy 
temperature Tp by means of evaporative cooling in an optical trap.^^^ The neutral atom of 
173Yb has sextuplet degeneracy owing to the degrees of freedom of nuclear spin I = 5/2 
with electron spins being quenched in the singlet state 5 = 0. Then, the sextet condensation 
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Fig. 13. (Color online) Temperature (T) dependences of bi{T) and b 2 {T) for fillings of fermions, 
^x/t = 0 and ^,/t = 3.9. 


is possible if a sufficiently attractive interaction works between two atoms in the dilute gas 
state. It has also been reported that the s-wave scattering length in the low energy limit of 
scattering atoms is positive and ~ 10.6 nm, which is fairly long compared with the range 
of a two-atomic interaction.This implies that there exists a shallow two-body bound state 
with the binding energy 


= -- 




mat 


(50) 


Then, according to the Nagaoka-Usui theorem,^®) the ground state of a six-particle system 
is fully symmetric in space coordinates and anti-symmetric in spin coordinates. This state 
is not an aggregation of two-atomic bound states, but is a coherent object formed by six 
particles. Of course, the situation is different in macroscopic systems.Nevertheless, there 
may be a chance that the sextet condensation is much more favorable than the Cooper pair 
condensation in some regions of temperature and atomic number density, as discussed in Sect. 


4. 


The binding energy, given by Eq. (50), with as = 10.6 nm, is estimated as |Eo|/A;b — 25 ^K. 
This is higher than the Fermi temperature Tp ~ 5 /rK of Yb gas attained from that with a 
temperature T ~ 100/iK and an atomic number density N/V ~ 7.3 x 10^^/cm^ at the initial 
stage of cooling. The Yb gas is finally cooled to T ~ 75 nK and N/V ~ 6.0 x 10^^/cm^ by 
evaporation. Thus, the cooling is accompanied by the dilution of the atomic number density, 
which decreases Tp. In the final stage of cooling, T ~ 0.37Tp, with Tp ~ 203nK. Therefore, 
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in the intermediate stage of cooling, there is a chance that both Tp and T of the system are 
comparable to or smaller than |£'o|/A:b- 

Here, let us estimate the strength of the attractive interaction potential, Vq, discussed in 
Sect. 4. We assume Vq as follows; 




(51) 


V, 0 < q < kc, 

0 , kc < q. 

Since Vq is a matrix element of scattering, (p, p') —)■ (p + q, p^ — q), and the scattering with 


IpI) IpI) Ip + hI) Ip^ — q| ~ fcp is important, it may be reasonable to take kc ~ kp. Then, the 


strength of the attractive interaction V* (< 0) in real space is related to V as 




1 ) 


(52) 


The strength of |H*| should be larger than the binding energy of the two-body bound 
state, i.e., \V*\ > |£'o|- Then, by using Eq. (52) and A^(eF) = mkp/27r‘^h‘^, 




(53) 


Therefore, it is really possible for the strong coupling region, |1 /|A^(£f) > 1, to be reached 
in the course of cooling and in the region of T where the sextet condensation is realized, as 
shown in Fig. 7. 

As discussed partly in Sect.4, the physical picture in the strong-coupling region is not sim¬ 
ple. The binding energy of the 6-body bound state is larger than that of three 2-body bound 
states, so that 6-atomic molecules are formed as T decreases. On the other hand, the Tbec 
of 6-atomic molecules is lower than that of diatomic molecules. Therefore, when the temper¬ 
ature is decreased from the normal state, the transition to the Bose-Einstein condensation of 
diatomic molecules would occur first if the diatomic molecules were formed at that tempera¬ 
ture. However, 6-atomic molecules are formed first when the temperature is decreased from the 
high-temperature side. Then, the formation of diatomic molecules is prohibited energetically, 
so that the Bose-Einstein condensation of diatomic molecules does not occur. 


8. Summary 

We have developed a mean-field theory for 2n-body (n > 2) condensation of the Ginzburg- 
Landau (GL) type, on the basis of the idea of variational principles on which the GL theory is 
based. We have found that the transition temperature Tc is expressed in concise form, which is 
numerically tractable for any number of n > 2. Namely, the Tc’s for the quartet, sextet, octet, 
and dectet condensations have been calculated for fermions with internal degrees of freedom. 
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4, 6, 8, and 10, respectively, not only in three-dimensional free space but also in a two- 
dimensional square lattice. We have also calculated the for the Cooper pair condensation 
with the same formalism of numerical calculations. The results are summarized as follows: 

1) There exists a threshold |Tth| of the strength of an attractive interaction V for the 
2n-body (n > 2) condensation to be realized, and the Tc’s exhibit the reentrant behavior for 
\V\ near the threshold |Vth|- In the case of three-dimensional free space, the threshold values 
extend as |Tth|-^(£F) = 0.1 ~ 0.3 from the dectet condensation to the quartet condensation. 
In the region of |Tth|-^(£F) in which 2n-body condensation has a finite T^, T^s are higher than 
that of the Cooper pair condensation. However, in the weak-coupling region |I4h|I\^(eF) < 0.1, 
the quartet condensation is not possible, while the Cooper pair condensation is always possible 
if V is attractive no matter how small |H|A^(eF) is. 

2) A similar trend is obtained in the case of a two-dimensional square lattice. A new 
aspect is the filling dependence of Tc for the quartet condensation. Tc increases as the filling 
UA increases. In the strong-coupling region U/t > 2.5, Tc’s are higher than that of the Cooper 
pair condensation for any filling 0 < ua < 1. The transition to the quartet condensed state 
is shown to be of the second order by an explicit calculation of the quartic terms in the GL 
thermodynamic potential. 

3) The sextet condensation is possible in a cold atom system of ^'^^Yb, which has a shallow 
two-body s-wave bound state implying that the attractive interaction satisfies the condition 
for the sextet condensation to occur dominating the Cooper pair condensation. 

Our GL-type formalism also makes it possible to search for the thermodynamic properties 
of 2n-body (n > 2) condensation near Tc, as discussed in Sect. 6.2 for the quartet (n = 2) 
condensation. However, detailed discussions for 2n-body (n > 3) condensation are left for 
future studies. The present results are valid near the transition temperature Tc because they 
are derived on the basis of the GL-type formalism. Therefore, it is not self-evident whether 
the 2n-body (n > 2) condensed state remains as the most stable ground state even though 
Tc’s are higher than that of the Gooper pairing state. Another important issue is how to treat 
the effect of the center-of-mass motion of 2n-body (n > 2) molecules in the strong-coupling 
regime, in order to discuss the crossover to the Bose-Einstein condensation of such molecules, 
as discussed by Nozieres and Schmitt-Rink in clarifying the problem of the BCS and BEC 
crossover phenomenon. 
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Appendix A: Ginzburg-Landau Formalism Revisited 

In this Appendix, we reformulate the Ginzburg-Landau (GL) theory^^^ for a uniform (s- 
wave spin singlet) pair condensed state by using the Feynman diagram representation. 

Let us start with the Feynman inequality for the thermodynamic potential 

n < flmf + (A'l) 

where H is the Hamiltonian of the system in consideration, is a mean-field Hamiltonian, 
and rimf is the thermodynamic potential for the system described by Rmf- Let us define the 
right-hand side of Eq. (A-1) as n, which is finally identified with the GL thermodynamic 
potential. Namely, 

H = riinf + {H — (A-2) 

The Hamiltonian of the fermion system with a pairing interaction Lk,k' is expressed as 
H = '^ CfcoLaka + X] X] k'®k+q/2,t®-k-|-q/2,i®-k'+q/2,;ak'+q/2,t> ("^'3) 

k,(T q k,k^ 

where is the dispersion of quasiparticles measured from the chemical potential, and 
(oko-) is the creation (annihilation) operator of quasiparticles with a wave vector k and a spin 
a (=t) i)- Hereafter, Vk^k' is assumed to be constant V {< 0). The mean-field Hamiltonian 
with the mean-field gaps Ak and is given by 

Lfmf = ^ ^ Cfc®ko-®k(T ~ ^ ^ ^^k®—k4,®kt T ^k^k'I'O'LkJ,^ • 

k,(j k 

Therefore, the operator corresponding to the second term in Eq. (A-2) is given by 

H - Rmf = ^ Lk,k'akt«-k.^“-k't«k't + X] (^k«-k.^akt + ^kakt«^k.^) • (^’5) 

k,k' k 

First, we calculate Hmf by perturbation expansion with respect to the s-wave gap A 
(without the k dependence) in the Hamiltonian (A-4) up to the quartic term in A and A*. 
The result is given by 

- Ho - Ai(r)|A|2+iA2(r)|A|4 + • • • . (A-6) 

Here, Hq is the thermodynamic potential in the normal state, the coefficients Ki{T) (i = 1, 2) 
are given by the Feynman diagrams shown in Fig. AT, and their analytical expressions are 
given as 

Ai(r) = T EE G(k,i6jG(-k,-i6„) 

en k 
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~ A^f J ^ tanh 


~ A^f log 


2ec7 


K 

2 


ttT 


and 


K2{T) =tY^ ien)G'(-k, -kO]' 

En k 

1 


tY,Nf 

Nf 7C(3) 


-e/^le+eir 


(A-7) 


(A- 8 ) 


(7rr)2 8 ’ 

where G is the Matsubara Green function of quasiparticles in the normal state, A^f is the 
density of states of quasiparticles at the Fermi level per spin, and £c is the energy cutoff of 
the pairing interaction. 7 is the Euler number and (^{z) is the Riemann (^-function. 



Fig. A-1. Feynman diagram for flmf up to the quartic terms in A and A*. 


Next, we calculate the grand canonical average of Eq. (A-5) with the mean-field Hamil¬ 
tonian (A-4) up to the quartic terms in the gaps A and A*. These terms are given by the 
Eeynman diagrams shown in Eigs. A2(a)~A2(e), and their analytical expressions are given as 

{H - ~ 2i^i(r)|A|2 + VKi{Tf\A\^-2VKi{T)K2{T)\A\^-2K2{T)\A\^ + ■■■ , (A-9) 


where the first term corresponds to Eig. A-2(a), the second term to Fig. A-2(b), the third term 
to Fig. A-2(c), and the fourth term to Fig. A-2(d), while the term in Fig. A-2(e) vanishes in the 
case where the particle-hole symmetry is maintained, as usually assumed in the weak-coupling 
treatment of the Cooper pair condensation. Indeed, an explicit expression for the triangle of 
the last term in Fig. A-2(e) is given as 


En k 


n 


1 1 

- ?k (-iCn - C-k)^ ’ 


(A-IO) 
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which is easily shown to be zero owing to the even-oddness of the integrand with respect to 
the inversion of e„ ^ —e„ and ^ —)• — 

Therefore, by adding Eq. (A-6), the GL thermodynamic potential 0(A) is expressed as 

0(A) ~ Oo + [1 + VK^{T)] [A:i(T)|A|2-2i^2(r)|A|4] + ii^2(T)|A|4 + • • • . (ATI) 

This is nothing but the GL thermodynamic potential. Indeed, 0(A) is exactly the same as 
Ogl(A) given by Leggett in Sect. 5.E of Ref. 27. The transition temperature T ^. is given by 
the condition that the coefficient of |Ap term is zero: 

1 = \V\Ki{T,). (AT2) 

By using Eq. (A-7), an explicit form of Eq. (AT2) is reduced to the BCS formula 

l = |R|iVFlog(^^). (ATS) 

Note that the quartic term of 0(A) is given essentially by the third term of Eq. (ATI) because 
the quartic term in the second term is not effective near the transition temperature Tc where 
the factor [1 -|- VKiiT)] vanishes. 



Fig. A-2. Feynman diagrams for {H — i7mf)mf up to the quartic terms with respect to A and A*. 
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The equation determining the gap A in the equilibrium at T < is given by the condition 
5n(A)/9A = 0, the explicit form of which is expressed as 

[1 + VKi{T)] Ki{Tc)A* + K2{Tc)\A\‘^A* = 0. (A-14) 

By using the explicit forms of Ki{T), i.e., Eq. (A-7), and that of K 2 {T), i.e., Eq. (A- 8 ), this 
equation is reduced to 

This is exactly the same form as that given by Gor’kov on the basis of the field theoretical 
method.^®’ 
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Appendix B: Calculation of G(r, r) 

In this Appendix, we calculate the quantity G(r,r) in the square brackets of Eq. (11); 

dp 


G(r,r) = 


i/(p)prE' 


1 


piP r 


(B-l) 


cn 

Considering the periodicity of the Matsubara Green function, we restrict the variable region 
of r within 0 < r < /3. Then, the summation with respect to is performed in a standard 
manner as 




1 


i^n Cp 


1 


dz 


27ri e ^2 + 1 z — 


e-<v 
g-/3?p _|_ 1 


Therefore, 


G(r,T) = 


do o ■ g(d 

^ |/(p)|Vp-(-i)" 


(27r)3^ ^'eKv + l 

After integrating with respect to the angular variables of p, we obtain 


G(r,r) = - 


1 1 

27r^ r 


dp |/(p)P p sin(pr) 


g(d-T)Sp 
g/3?p -P 1 


(B-2) 


(B-3) 


(B-4) 


Changing the integration variable from p to x = = {p^/2meY) — 1, and using the form 

of /(p) given by Eq. (25), we obtain 


G(r,r) = - 


m Ep 
27r^ 


r 


/Xc 

dx sin[\/WT(A:irr)] 


g(/ 3 -T)£Fa: 
ePepx I ’ 


(B-5) 


where Xc = (ec/ep) — I, and we have used an approximation /x ~ ep. 
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Appendix C: Calculations of quartic terms in A and A* 

In this Appendix, we give the explicit expressions of the quartic terms in A and A* of 
the GL expansion for thr quartet condensation in a two-dimensional square lattice where the 
function /(pi) is set to unity, i.e., /(pi) = 1 . 


rimf CLS function of A and A* 


Of the quartic terms in A and A*, those for are given by the Feynman diagrams 
shown in Fig. C-1. 

The analytical expression C'i(r) for the diagram shown in Fig. CT(a) is given as 

Ci(r) = n T ^ G(p., i£„.) 

i=l Pi Cnj 

X(5 (pi + P2 + P3 + P4) X 
x6 (pi + P 2 + P5 + Pe) X 

X(5 (P5 P6 P7 + Ps) X (C-1) 

Here, the combination factor — 4 C 2/2 comes from the number of combinations for perturbation 
expansion and the Wick theorem. Namely, 

X 4 C 2 X 2 X 4 C 2 X (-hi) =-^, (C-2) 

where the factor (—1/4!) comes from the perturbation expansion of Hmf to the 4th order in A 
and A*, the factor 4 C 2 is the number of ways of choosing two A’s from four products of the 
perturbation terms in H^f, given by Eq. (3), the factor 2 is the number of ways of choosing two 
A*’s from two products of the perturbation terms in also given by Eq. (3), another factor 
4 C 2 is the number of combinations how to choose 2 spin states of Green functions connecting a 
certain pair of A and A* from a, j 3 , 7 , and 6 and the factor (-1-1) represents that the number 
of interchanges of Eermion operators is even in the Wick expansion. Other assignments of 
the spin variables a, ft, 7 , and <5 to Green functions are automatically determined by the 
conservation law of spins. 

In the case of lattice systems, instead of the relation Eq. (9), the following relation holds: 

(Pl + P2 + P3 + P4) = ^ V e'(Pl+P^+P3+P4)-r.^ (C-3) 

X 

where Nl is the number of lattice points. Then, by using Eqs. (G-3) and (10), the coefficient 
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2=1 ' 


n 2 


Nl 


^r^G(qi,iemi)e' 


(qi-ri —emiTi) 


Q.1 £77 




J{^2■r2-em2'r2) 


<12 




H 2 


+r3)-em3 (n+rs)) 


iVu 


qs fmo 


X!/ ^ X!/ <^(^4, iem4) 


„i(q4'(r2+r3) ^m4 (T2+r3)) 


(C-4) 


<l4 

This expression is managed easily using the FFT algorithm, as discussed in Sect. 5. In 
terms of Xm defined in Eq. (38), (C-4) is expressed as 

Ci{T) = ^'^ 2 /' dr3X2(rf\ri)^2(rf\T2) 

^ 

^ i ^ i ^ i 


xX 2 (rf^ + rfVi + 'r3)X2(rf^ + rf\r 2 + rs). (C-5) 

Then, by calculations similar to those leading to Eq. (44) from Eq. (42), the expression (C-5) 
is reduced to 


Ci{T) = j;[X2(-k,-icu„)X2(k,ia;„)]2, (C- 6 ) 

where ujn is the bosonic Matsubara frequency. Note that Xm(rj,r) with the odd natural 
number m is expanded into the Fourier series with the component Xm(ri, ie^) with a fermionic 
Matsubara frequency €n = (2n + l) 7 rr because Xm{ri,T + P) = —Xm(ri,r) for the odd natural 
number m. 

The analytical expression C 2 {T) for the Feynman diagram shown in Figs. C-l(b) is given 


C 2 (r) = n E G(p., i£„,) 

i=l Pi En; 

x5 (pi + P 2 + P 3 + P 4 ) ^ + 6712 +£ 713+£714 ,0 

x 5 (p 2 + P 3 +P 4 + P 5 ) ^ ^^ 712 +^773+€77.4+^715,0 

x5(p5 + P6+P7 + P 8 ) ^ ^€7i3+e7T,g+eTT,y+e7ig,0* (C-7) 


Here, the combination factor 4 C 3/2 comes from the number of combinations for perturbation 
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expansion and the Wick theorem. Namely, 

X 4 C 2 X 2 X 4 C 3 X (-1) = (C- 8 ) 

where the factor (—1/4!) comes from the perturbation expansion of flmf to the 4th order, the 
factor 4 C 2 is the number of ways of choosing two A’s from four products of the perturbation 
terms in given by Eq. (3), the factor 2 is the number of ways of choosing two A*’s 

from two products of the perturbation terms in also given by Eq. (3), the factor 4 C 3 

is the number of combinations how to choose 3 spin states of Green functions connecting a 
certain pair of A and A* from a, /?, 7 , and 5, and the factor (- 1 ) represents that the number of 
interchanges of Fermion operators is odd in the Wick expansion. Other assignments of the spin 
variables a, j3, 7 , and 5 to Green functions are automatically determined by the conservation 


law of spins. 

By using Eqs. (C-3) and (10) and similar ones, the coefficient C 2 {T) is reduced to 


JL rP 




„i(qi-ri-e™^ri) 




r2-<imoT2) 


Q2 ^^^2 


qs Cmg 

Q4 


In terms of Xm defined in Eq. (38), the right-hand side of Eq. (C-9) is expressed as 
C'2(T) = ^^^^/ dn f dT2 I dr3Xi(rfVi)^3(rS^\T-2) 


,.( 1 ) ^( 2 ) „( 3 ) 

^ i i i 


/o Jo Jo 


xA 3 (rf^ +rf\T-i +'r 3 )-Ei(rf^ -krf\r 2 + T 3 ). 


(G-10) 


Then, by calculations similar to those leading to Eq. (44) from Eq. (42), the expression (GTO) 
is reduced to 

C2(r) = E E ’ (c-ii) 

k en 

where is the fermionic Matsubara frequency as mentioned just below Eq. (C- 6 ). 

The analytical expression C^{T) for the Feynman diagram shown in Fig. GT(c) is identical 
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Cs{T) = C2{T). 


(C-12) 



Fig. C-1. Feynman diagrams for rimf of the quartic terms with respect to A and A*. 


Terms without V in {H — 

Of the quartic terms in A and A* arising from {H — those without the interaction 

V are given by the Feynman diagrams shown in Fig. C-2. 

The analytical expression Di(T) for the diagram shown in Fig. C-2(a) is given as 

8 1 

Di{T) = 2XiC2T^ G(p*,ie„J 

Pi 

Xd (pi + P2 + P3 + P4) X 
Xd (pi + P2 + P5 + Pe) X 

Xd (p5 + P6 + P7 + Ps) X (C-13) 

This is the same as Eq. (CT) except for a difference in a prefactor. Here, the factor 2 x 4 C 2 
comes from the number of combinations for perturbation expansion and the Wick theorem. 
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Namely, 

^x2x3x2x 4 C 2 X (+1) = 2 X 4 C 2 , (C-14) 

where the factor (1/3!) comes from the perturbation expansion of {H — -ffmf) to the 3rd order 
in A and A*, the factor 2 is the number of ways of choosing A or A* from [H — 
given by Eq. (5), the factor 3 is the number of ways of choosing two A*’s or A’s from three 
products of the perturbation terms in given by Eq. (3), another factor 2 is the number 
of ways of choosing two A*’s from two products of the perturbation terms in also given 
by Eq. (3), the factor 4(72 is the number of combinations how to choose 2 spin states of Green 
functions connecting a certain pair of A and A* from a, I3, 7 , and <5, and the factor (+1) 
represents that the number of interchanges of Eermion operators is even in the Wick expansion. 
Other assignments of the spin variables a, I3, 7 , and S to Green functions are automatically 
determined by the conservation law of spins. Therefore, Di{T) is given in terms of C'i(T) as 

Di{T) = -4Ci(T). (G-15) 

The analytical expression D2{T) for the Eeynman diagram shown in Eig. GT(b) is given 
as 

® 1 

D2{T) = -2x,C,T^ll — J212 G(Po i 6 nj 

i=l Pi 

x5 (pi + P2 + P3 + P4) X 
x5 (P2 + P3 + P4 + Ps) X 

x5 (P5 + P6 + P7 + Ps) X (C-16) 

This is the same as Eq. (G-7) except for a difference in a prefactor. Here, the factor —2 x 4 G 2 
comes from the number of combinations for perturbation expansion and the Wick theorem. 
Namely, 

^x 2 x 3 x 2 x 4 G 3 X (-1) = -2 X 4 G 3 , (G-17) 

where the factor (1/3!) comes from the perturbation expansion of {H — to the 3rd order 
in A and A*, the factor 2 is the number of ways of choosing A or A* from [H — 77mf)) 
given by Eq. (5), the factor 3 is the number of ways of choosing two A*’s or A’s from three 
products of the perturbation terms in i7mf; given by Eq. (3), another factor 2 is the number 
of ways of choosing two A*’s from two products of the perturbation terms in i7mf; also given 
by Eq. (3), the factor 4(73 is the number of combinations how to choose 3 spin states of 
Green functions connecting a certain pair of A and A* from a, /3, 7 , and 5, and the factor (- 1 ) 
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represents that the number of interchanges of Fermion operators is odd in the Wick expansion. 
Other assignments of the spin variables a, /3, 7 , and S to Green functions are automatically 
determined by the conservation law of spins. Therefore, D 2 {T) is given in terms of 0 * 2 (T) as 

D2{T) = -AC2{T). (C-18) 

The analytical expression D^{T) for the Feynman diagram shown in Fig. C-2(c) is identical 
to that shown in Fig. C-2(b). Namely, 

Ds{T) = D2{T). (C-19) 

Equations (CT5), (CT 8 ), and (CT9) indicate that the contributions of Fig. C-2 are twofold 
those of Fig. CT in size and opposite in sign. 



Fig. C-2. Feynman diagrams for (i/— of the quartic terms with respect to A and A* without 

the interaction V. 


Terms with V in {H — 

The quartic terms including the interaction V in {H — are given by the Feynman 

diagrams shown in Fig. C-3. 
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The analytical expression Di{T) for the diagram shown in Fig. C-3(a) is given as 

i —1 Pi 

x5 (pi + P2 + P7 + Ps) X 

X(5 (p7 + P8 - P9 - Pio) X 
X(5 (P5 + P6 + P9 + Pio) X 

X(5 (P3 + P4 + P5 + Pe) X (C-20) 


Here, the combination factor 4 C 2 comes from the number of combinations for perturbation 
expansion and the Wick theorem. Namely, 


4! 


X 4 C 2 X 4 C 2 X 2 X 2 X (+1) = 4 C 2 , 


(C-21) 


where the factor (1/4!) comes from the perturbation expansion of the first term of Eq. (5) to 
the 4th order in A and A*, the factor 4 C 2 is the number of ways of choosing 2 spin states in 
the interaction V from a, 13, 7 , and 6, another factor 4 C 2 is the number of ways of choosing 
two A’s from four products of the perturbation terms in Fimf, given by Eq. (3), the factor 
2 X 2 is a product of the number of ways of choosing two A’s from two products of the 
perturbation terms in H^f, also given by Eq. (3), and that of choosing A*, and the factor 
(+1) represents that the number of interchanges of Eermion operators is even in the Wick 
expansion. Other assignments of the spin variables a, (3, 7 , and <5 to Green functions are 
automatically determined by the conservation law of spins. 

By using Eqs. (C-3) and (10) and similar ones, the coefficient D^iT) is reduced to 

r n 2 


Da{T) 

V 



^ qi Cmi 


H 2 


Nl 




<12 ^^712 




7(q3'(-r2+r3)-e™3(-T2+T3)) 


Nl 


qs fn 


q 2 


E^E G'(q4,ie^4)e'('i"'("3+"")-""‘4 03+r4)) 


<l4 ^m4 


E ^ E ^(qs, 


<l5 ^^5 


(C-22) 
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By calculations similar to those leading to Eq. (C-6) [Eq. (C-11)] from Eq. (C-4) [Eq. (C-9)], 
the expression (C-26) is reduced to 

EE [X2(-k,-iw„)X2(k,ia;n)]^ (C-23) 

k LOn 

where ujn is the bosonic Matsubara frequency. 

The analytical expression D^{T) for the diagram shown in Fig. C-3(b) is given as 

^ = -2 X 4Cj U n T y; ^ G(p., i,,._) 

i —1 Pi 

X5 (pi + P5 + P6 + pg) X 
X5 (p5 + P6 - P7 - Ps) X 
X5 (p7 + P8 + P9 + Pio) X 

X5 (p 2 + P3 + P4 + Pio) X (C-24) 

Here, the combination factor — 2 X 4 C 2 comes from the number of combinations for perturbation 
expansion and the Wick theorem. Namely, 

^ X 4 C 2 X 4 C 2 X 2 X 2 X (—1) = —2 X 4 C 2 , (C-25) 

where the factor (1/4!) comes from the perturbation expansion of the first term of Eq. (5) to 
the 4th order in A and A*, the factor 4 C 2 is the number of ways of choosing 2 spin states in 
the interaction V from a, 13, 7 , and 6, another factor 4 C 2 is the number of ways of choosing 
two A’s from four products of the perturbation terms in given by Eq. (3), the factor 2x2 
is a product of the number of ways of choosing two A’s from two products of the perturbation 
terms in also given by Eq. (3), and that of choosing A*, the factor 2 is the number of 
ways of choosing a spin state, S or 7 , for the Green fnnction on the left side of Fig. C-3(b), 
and the factor (-1) represents that the number of interchanges of Fermion operators is odd in 
the Wick expansion. Other assignments of the spin variables a, f3, 7 , and 6 to Green functions 
are automatically determined by the conservation law of spins. 

By using Eqs. (G-3) and (10), the coefficient D^(T) is reduced to 

= -2XiCiH [’’‘‘nY. /-E’’EG'(qi.iE™.)e'">‘-''‘— 

i=\ ^ r? Ql €.mi 

X ^ E ^ E G'(q2, 

<l2 ^m2 
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H 2 


Nl 




(ri+r4)-e^3 (Ti+r4)) 


qs c„,3 




Nl 


1714 ) 


^i(q4-(r2-r4)-em4 (^- 2 -^ 4 )) 


0,4 ^rriA 




05 £7715 


(r 2 +r 3 )-€m 5 (T 2 +T 3 )) 


n 3 


I • (C-26) 

06 ^Tng 

By calculations similar to those leading to Eq. (C- 6 ) [Eq. (C-11)] from Eq. (C-4) [Eq. (C-9)], 
the expression (C-26) is reduced to 

^5(^) ^ 

2x 4C2^ E E [-^i(-ki,-ieni)^2(k2,i‘^n2)]^-’^i(ki-k2,ie„3-ia;n2)^3(ki,ieni), 

^ ki,k2 €n^ ,^712 


V 


(C-27) 

where ujn and are the bosonic and fermionic Matsubara frequencies, as mentioned just 
below Eq. (C- 6 ). 

The analytical expression Dq{T) for the diagram shown in Fig. C-3(c) is given as 


mT) 

V 


= -4C2 


n jVi ^ ^ ^ 

i=l Pi e-n. 


X(5 (pi + P2 + P3 + P 4 ) X (5, 


^ni H“^n2 H”^n4 jO 


(C-28) 


Here, the combination factor — 4 C 2 comes from the number of combinations for perturbation 
expansion and the Wick theorem. Namely, 

1 


4! 


X 4 C 2 X 4 C 2 X 4 X 2 X 2 X (—1) = — 4 C 2 , 


(C-29) 


where the factor (1/4!) comes from the perturbation expansion of the hrst term of Eq. (5) to 
the 4th order in A and A*, the factor 4 C 2 is the number of ways of choosing 2 spin states in 
the interaction V from a, 13, 7 , and 6, another factor 4 C 2 is the number of ways of choosing 
two A’s from four products of the perturbation terms in given by Eq. (3), the factor 
2 X 2 is a product of the number of ways of choosing two A’s from two products of the 
perturbation terms in Hy^i, also given by Eq. (3), and that of choosing A*, and the factor (-1) 
represents that the number of interchanges of Fermion operators is odd in the Wick expansion. 
Other assignments of the spin variables a, 13, 7 , and 5 to Green functions are automatically 
determined by the conservation law of spins. 
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Let us define the quantity in the square brackets in Eq. (C-28) by Z(T). By using Eqs. 
(C-3) and (10), Z(T) is reduced to 


rP 

Z{T) = T^ drY" 

Jo ^ 




Ql £mi 


n 2 


1 


(C-30) 


<12 ^rn2 

By calculations similar to those leading to Eq. (C-6) [Eq. (C-11)] from Eq. (C-4) [Eq. (C-9)], 
the expression (C-30) is reduced to 


T 


Xsi-k, -ie„)[Xi(k,ie„)]' 


(C-31) 


k en 


Then, Dq{T), given by Eq. (C-28), is given as 

DoiT) 


V 


= -4C2 [Z{T)Y 


(C-32) 


The analytical expression Dj{T) for the diagram shown in Eig. C-3(d) is given as 

10 


DijT ) 
V 




■ 1 

1=1 Pi En; 

X(5(pi + P7 + P9 + Pi0) ^ +en7+en9+€njQ ,0 


X(5 (P2 + P8 - P5 - P?) X (5, 




X(5(p6 + P8+P9 + Plo) X +eng+en9+€n]^Q ,0 
X(5 (p 3 + P4 + P5 + P6) ^ ^€^3+en4+€n5+€ng ,0- 


(C-33) 


Here, the combination factor 4 C 2 comes from the number of combinations for perturbation 
expansion and the Wick theorem. Namely, 

^ X4 C 2 X4 C 2 X 2 X 2 X (+1) = 4 C 2 , (C-34) 

where the factor (1/4!) comes from the perturbation expansion of the first term of Eq. (5) to 
the 4th order in A and A*, the factor 4 C 2 is the number of ways of choosing 2 spin states in 
the interaction V from a, 13, 7 , and 6, another factor 4 C 2 is the number of ways of choosing 
two A’s from four products of the perturbation terms in H^ai, given by Eq. (3), the factor 
2 X 2 is a product of the number of ways of choosing two A’s from two products of the 
perturbation terms in Hmi, also given by Eq. (3), and that of choosing A*, and the factor 
(+1) represents that the number of interchanges of Eermion operators is even in the Wick 
expansion. Other assignments of the spin variables a, j3, 7 , and 6 to Green functions are 
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Fig. C-3. Feynman diagrams for {H — of the quartic terms with respect to A and A*, which 

include the interaction V. 


automatically determined by the conservation law of spins. 

By using Eqs. (C-3) and (10), the coefficient D^iT) is reduced to 

4 


Di{T) 

V 


= 4C2 n / dn Y, 

i=l Ti 


2 




Nl 


0,1 Cmi 


Nl 


Nl 




+r2)-em2 (ri+'r2)) 


O 2 ^rn2 


E^E G(q3,ie^3)e^('l3-(ri+r4)-.^3(ri+r4)) 


03 ^^3 
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Nl 


J]Tj^G'(q4,ie„ 


„i{q4-(r2-r4)-em4 (T2-r4)) 


Q4 ^m 4 


H 2 




Q5 ^^5 


^ X!/ ^(^ 6 , iCme; 


oi{q5'(r2+r3) ^mc (T 2 +T 3 )) 


oUqe-rs-emg'Ts) 




<16 ^^6 


E^E 


QT €7717 


=i(q8'(-r4)-emo(-r4)) 


(C-35) 


<18 ^'^8 

By calculations similar to those leading to Eq. (C- 6 ) [Eq. (C-11)] from Eq. (C-4) [Eq. (C-9)], 
the expression (C-35) is reduced to 


MT) 


rjn4 


y =4C2^X] E E ^i(-ki,-ieni)^i(ki-k3,ieni-iwng) 

ki~k 4 ^ni }^ n2 ^n 3 5^714 

xXi(-ki + ka - k4, -kni + iWns “ iWn 4 )Xi(-k 2 , -iCnj) 

X Xi (k2 - ks, i€n2 - i‘^n3 ) -^2 (ks, iw^g). 


xXi(-k 2 + ka - k4, -kna +iWn3 - ia;n4)X2(k4,iWn4)- (C-36) 

The analytical expression D^{T) for the diagram shown in Fig. C-3(e) is given as 

^ = 2x4c,un^EE0(p..i^».) 

i—1 Pi 

X(5 (pi + P7 + P8 + Pio) X 
X(5 (P7 + P8 - P2 - Pg) X 
X(5 (P5 + P6 + P9 + Pio) X 

x-5 (P3 + P4 + P5 + Pe) X (C-37) 

Here, the combination factor 2 x 4 C 2 comes from the number of combinations for perturbation 
expansion and the Wick theorem. Namely, 

^ X 4 C 2 X 4 C 2 X 2 X 2 X 2 X (+1) = 2 X 4 C 2 , (C-38) 

where the factor (1/4!) comes from the perturbation expansion of the first term of Eq. (5) to 
the 4th order in A and A*, the factor 4 C 2 is the number of ways of choosing 2 spin states in 
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the interaction V from a, (3, 7 , and 5, another factor 4 C 2 is the number of ways of choosing two 
A’s from the four products of the perturbation terms in given by Eq. (3), the factor 2x2 
is a product of the number of ways of choosing two A’s from two products of the perturbation 
terms in H^f, also given by Eq. (3), and that of choosing A*, the factor 2 is the number of 
ways of choosing a spin state, S or 7 , for the Green function on the left side of Fig. C-3(e), and 
the factor (+1) represents that the number of interchanges of Fermion operators is even in 
the Wick expansion. Other assignments of the spin variables a, 13, 7 , and 6 to Green functions 
are automatically determined by the conservation law of spins. 

By using Eqs. (C-3) and (10), the coefficient Ds{T) is reduced to 


DsjT) 

V 


/9 


2X4C2 n / 


2=1 ' 


1 






Ql €.mi 

i(q2-(ri+r2)-em2 ('ri+'f2)) 


<12 


- 2 


Nl 


<13 ^mq 




<14 £7714 




(r2-r4)-em^ (r2-T4)) 


(r2+r3)-em5 (T2+r3)) 


<15 ^^5 




n 2 


rs-£ 7776 ^- 3 ) 


Nl 


<16 ^^6 




„i(q6-(-r4)-em6(-T-4)) 


(C-39) 


<l 7 £7777 

By calculations similar to those leading to Eq. (G- 6 ) [Eq. (G-11)] from Eq. (C-4) [Eq. (G-9)], 
the expression (G-39) is reduced to 


V 


^ — 2 X 4G2-^ ^ ^ ^ Ai(ki,ie„j)Xi(—ki, —ienj^i(ki — k2,ie„ 

^ ki~k 3 £77]^ ^772 5^773 

xXi(-ki + k2 - ks, -iCni + - iu;n3)X2(k2,iuJn2) 

xX2(k3, ia;„3)A2(—k3, —ia;„3). 


(G-40) 
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The analytical expression Dg{T) for the diagram shown in Fig. C-3(f) is given as 

i=l Pi 

X5 (pi + P 2 + P3 + P4) X (5,„^+e„2+€„3+€„4,0 

X(5 (p 8 + P 9 - P 3 - Pio) X 
X(5 (pi + P 2 + P 8 + pg) X 

X(5(p4 + P5+P6 + P7) ^ ^€ 77,4 +£715 +€ng+€ 71 . 7 , 0 * (C-41) 

Here, the combination factor 4 C 2 comes from the number of combinations for perturbation 
expansion and the Wick theorem. Namely, 

^ X 4 C 2 X 4 C 2 X 2 X 2 X (+1) = 4 C 2 , (C-42) 

where the factor (1/4!) comes from the perturbation expansion of the first term of Eq. (5) to 
the 4th order in A and A*, the factor 4 C 2 is the number of ways of choosing 2 spin states in 
the interaction V from a, 13, 7 , and 6, another factor 4 C 2 is the number of ways of choosing 
two A’s from four products of the perturbation terms in Fimf, given by Eq. (3), the factor 
2 X 2 is a product of the number of ways of choosing two A’s from two products of the 
perturbation terms in H^f, also given by Eq. (3), and that of choosing A*, and the factor 
(+1) represents that the number of interchanges of Eermion operators is even in the Wick 
expansion. Other assignments of the spin variables a, (3, 7 , and <5 to Green functions are 
automatically determined by the conservation law of spins. 

By using Eqs. (C-3) and (10), the coefficient DgiT) is reduced to 

^ = 40 + fdnY: 

i=l ri ^ qi 

Q 2 

G'(q 3 , 02-^4)) 

<13 

X -^Y^Y G'(q4, 

<14 €m4 

<15 ^^5 
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^ . (C-43) 

Q 6 ^^6 

By calculations similar to those leading to Eq. (C-6) [Eq. (C-11)] from Eq. (C-4) [Eq. (C-9)], 
the expression (C-43) is reduced to 

'^^=4^2^E E [^i(-ki,-ienJ]'X2(k2,ia;„,)X2(-k2,-io;„J 

^ kl ,k2 :^n2 

xXi(ki -k2,ieni - ia;„2)X3(ki,ieni)- (C-44) 
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